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Skyrmions are topological spin textures of interest for fundamental science and applications. 
Previous theoretical studies have focused on systems with broken bulk inversion symmetry, where 
skyrmions are stabilized by easy-axis anisotropy. We investigate here systems that break surface- 
inversion symmetry, in addition to possible broken bulk inversion. This leads to two distinct 
Dzyaloshinskii-Moriya (DM) terms with strengths D±, arising from Rashba spin-orbit coupling 
(SOC), and Dy from Dresselhaus SOC. We show that skyrmions become progressively more sta¬ 
ble with increasing D±/D\\, extending into the regime of easy-plane anisotropy. We find that the 
spin texture and topological charge density of skyrmions develops nontrivial spatial structure, with 
quantized topological charge in a unit cell given by a Chern number. Our results give a design prin¬ 
ciple for tuning Rashba SOC and magnetic anisotropy to stabilize skyrmions in thin films, surfaces, 
interfaces and bulk magnetic materials that break mirror symmetry. 


Recently there has been a surge of interest in skyrmions 
in chiral magnetic materialJ^^^^, ranging from fundamen¬ 
tal science to potential device applications. A skyrmion 
is a spin texture characterized by a topological invariant 
that, in metallic magnets, gives rise to the topological 
Hall effect!^ and may also have implications for non- 
Fermi liquid behavior^. The ability to write and erase in¬ 
dividual skyrmion^, along with their topological stabil¬ 
ity, small size, and low depinning current densit}!^, paves 
the way for potential information storage and processing 
applications. 

Experiments have focussed primarily on skyrmions in 
non-centrosymmetric crystals with broken bulk inversion 
symmetry, e.g., metals like MnSi, FeGe and insulators 
like Cu 20 Se 03 . In these bulk materials, the skyrmion 
crystal (SkX) phase is stable only in a very limited re¬ 
gion of the magnetic field (77), temperature (T) phase 
diagrauPK^. On the other hand, the skyrmion phase is 
found to be stable over a much wider region of (T, 77) in 
thin films of the same materi al even extending 
down to T = 0 in some casefl^^. (A class of two di¬ 
mensional (2D) systems shows atomic-scale skyrmion^I^I 
arising from competing local interactions, distinct from 
the spin-orbit induced chiral interactions that we focus 
on here.) 

A key question that we address is this paper is: How 
can we enhance the domain of stability of skyrmion spin 
textures? We are motivated in part by the thin film ex¬ 
periments, and also by the possibility of chiral magnetism 
in new 2D systems like oxide interface^^^^^. We show 
how the SkX become progressively more stable over ever 
larger regions in parameter space of field 77 and mag¬ 
netic anisotropy A, as the effects of broken surface inver¬ 
sion dominate over those of broken bulk inversion. The 
key parameter responsible for this behavior is the ratio 
77±/77|| of the strength of the chiral magnetic interaction 
arising from broken bulk inversion (77||) to that arising 
from broken surface inversion (77^); see Fig. 1. 

We begin by summarizing our main results, which re¬ 
quires us to introduce some terminology. We focus on 


magnets in which spin textures arise from the inter¬ 
play between ferromagnetic exchange J and the chiral 
Dzyaloshinskii-Moriya (DM) interaction • (S^ x Sj). 
Spin-orbit coupling (SOC) determines the magnitude of 
the D vector, while symmetry dictates its direction. Bro¬ 
ken bulk inversion symmetry (r—r) leads to the Dres¬ 
selhaus DM term with = 77|| where Vij = | 

with Tij = (r^ —Tj). On the other hand, broken sur¬ 
face inversion or mirror symmetry (z^—z) leads to the 
Rashba DM term with Dij = D±(z x r^j). In the limit 
of weak SOC, B/J < 1, where 77 = (77^ + 77 ^)V2^ the 
length scale of spin textures is {J/D)a ^ a (the micro¬ 
scopic lattice spacing) and we can work with a continuum 
“Cinzburg-Landau” field theory. 

We show in Fig. 1 the evolution of the T = 0 phase di¬ 
agram going from the pure Dresselhaus limit to the pure 
Rashba limit. Each phase diagram is plotted as a func¬ 
tion of the (dimensionless) field 77J/77^ and anisotropy 
AJ/77^. Here A > 0 (A < 0) corresponds to easy-plane 
(easy-axis) anisotropy. Our main results are: 

(1) As the Rashba 77^ is increased relative to the Dres¬ 
selhaus 77||, the spiral and skyrmion phases become in¬ 
creasingly more stable relative to the vertical cone phase, 
and penetrate into the easy-plane anisotropy side of the 
phase diagram. 

(2) With increasing 77x/77||, the textures change con¬ 
tinuously from a Bloch-like spiral to a Neel-like spiral. 
Correspondingly, the skyrmion helicity evolves with a 
vortex-like structure in the Dresselhaus limit to a hedge¬ 
hog in the Rashba limit, which is shown to impact the 
ferrotoroidic moment. 

(3) In the pure Rashba limit, we find the largest do¬ 
main of stability for the hexagonal skyrmion crystal. In 
addition we also find a small sliver of stability for a square 
skyrmion lattice, together with an elliptic cone phase, 
distinct from the well-known vertical cone phase in the 
Dresselhaus limit. 

(4) We see in Fig. 2 that in the Rashba limit the spin 
texture of the skyrmion and their topological charge den- 
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FIG. 1. Phase diagrams as a function of AJ/D^ and HJ/D^ for four values of D±/D\\. Easy-axis anisotropy corresponds to 
A < 0 while easy-plane to A > 0. The cone, elliptic cone, and tilted FM phases are shown schematically, with the Q-vector 
shown in red and the texture traced out by spins shown in black. The color bar on the right indicates ruz for the elliptic cone 
and tilted FM phases in the = 0 panel. Insets: Unit cell in the hexagonal (Hex) skyrmion crystal (SkX) phase with white 
arrows indicating the projection of magnetization on the x-y plane. The colors indicates the magnitude and direction of the 
spin projection following the convention of ref. [3] indicated in the color wheel. Thick lines denote continuous transitions, while 
thin lines indicate first-order phase transitions. 


sity x(r) begins to show non-trivial spatial variations as 
one changes anisotropy, but the total topological charge 
^sk = f G^^rx(r) in each unit cell remains quantized, 
even when x(i*) seems to “fractionalize” with positive 
and negative contributions within a unit cell. 

(5) For H > 2M, one can have isolated skyrmions in 
a ferromagnetic (FM) background, and their topological 
charge A^sk is quantized, as usual, by the homotopy group 
= Z. For H < 2M, we find that skyrmions cannot 
exist as isolated objects, and A^sk must now be defined by 
the Z Chern number classifying maps from the SkX unit 
cell, a two-torus to 5'^, the unit sphere in spin-space, 
a definition that works for all values of Hj^A. 

Free energy: We consider a continuum (free) energy 
functional F[m] = f d^rJ^(m) with 

T = Tj ^ Jdm + — Hruz (1) 

whose form is dictated by symmetry. The isotropic ex¬ 
change term Tj = (^/2) 5])^(VmQ,)^ {a = x,y,z) con¬ 
trols the gradient energy through stiffness J . The DM 
contribution in the continuum 

= D cos f3 m-(Vxm)+D sin/3 m-[(zx V) xm] (2) 

is the sum of the two terms discussed above. The 
= DcosP term arises from Dresselhaus SOC in the 
absence of bulk inversion and = D sin p from Rashba 
SOC with broken surface inversion. The anisotropy term 
Ta = Arnl can be either easy-axis {A < 0) or easy-plane 
{A > 0). Several different mechanisms contribute to A, 
including single-ion and dipolar shape anisotropies. In 
addition, Rashba SOC naturally leads to an easy-plane, 
compass anisotropy A± D\l2J ^ which is energetically 


comparable to the DM terrrP^^. We treat M as a free, 
phenomenological parameter. 

We focus on T = 0 where the local magnetization is 
constrained to have a fixed length m^(r) = 1, and it 
should be hardest to stabilize skyrmions; once |m(r)| can 
become smaller due to thermal fluctuations, skyrmions 
should be easier to stabilize. It is convenient to scale all 
distances by the natural length scale J jD (setting the mi¬ 
croscopic a = 1) and scale the energy T by j J. All our 

results will be presented in terms the three dimension¬ 
less parameters that describe namely field HJ 
anisotropy AJ^ and tan^^ = 

Phase Diagram: In Fig. 1, we show the evolution 
of the (A^H) phase diagram as a function of tan^^ = 
Dx/D||, increasing from left to right. These results were 
obtained by minimizing the energy functional *Fdm sub¬ 
ject to m^(r) = 1. The energies of the fully polarized 
ferromagnet (FM), the tilted FM, and the vertical cone 
states can be determined analytically, while the energies 
for the spiral, the skyrmion crystals and the elliptic cone 
state were found by a numerical, conjugate gradient min¬ 
imization approach. In all cases, the results were checked 
by semi-analytical variational calculations. Details of the 
methodology are described in the Supplementary Mate¬ 
rials; here we focus on the results. 

We begin with well knowrP^ results in the Dressel¬ 
haus limit (left panel of Fig. 1), where the hexagonal 
SkX and spiral phases are stable only in a small region 
with easy-axis anisotropy (A < 0). The A > 0 region is 
dominated by the vertical cone phase, where mcone(^) = 
(cos(/?( 2 :) smOQ,smip{z) sincos ^o) with (p{z) = D\\z/ J 
and cos6>o = H/[2A + Dy/J]. The phase boundary be- 
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FIG. 2. Evolution of the spin texture m (top row) and the topological charge density x (bottom row) for four values of AJ/ 
at fixed HJ/D^ = 0.7 in the Rashba limit = 0). White arrows indicate the projection of m into the x-y plane. The colors 
also indicates the magnitude and direction of the spin projection following the convention of ref. [3] indicated in the color wheel. 
The development of nontrivial spatial variation in x(r) is discussed in the text. Note, however, that in each case integral over 
a single unit cell f d^rx(r) = — 1. 


tween the vertical cone and polarized FM is given by 

H = 2A^Dl/J. 

We note a change of variables that greatly simplifies 
the analysis of skyrmion crystal and spiral phases. This 
transformation is useful when m = m{x,y) has no z- 
dependence (along the field). Using a rotation Rz{—f3) 
by an angle —P about the z-axis, we define rv{x^y) = 
Rz{—P)m{x^y). It is then easy to show that the two 
terms in Q combine to give a pure Dresselhaus form 
Rbm = Dn {\/ X n). The other terms in 0 are invariant 
under this transformation, and thus R greatly simplifies 
using the transformation m^n. 

We choose, without loss of generality, x as the prop¬ 
agation direction for the spiral of period L, so that 
nsp(^) = ( 0 , sin 6 >(x), cos 6 >(x)) with nsp(x + I/) = nsp(x). 
(Note that this is not in general a single-q spiral). We 
minimize T to find the optimal L and optimal func¬ 
tion 0{x), which is a ID minimization problem. For 
the SkX, we first pick a unit cell, hexagonal or square. 
We then find its optimal size and optimal texture n = 
(cos cp sin sin (p sin cos 0 ), by solving a 2D minimiza¬ 
tion problem to determine (^(x, y) and 0{x, y) within a 
unit cell, subject to periodic boundary conditions. We 
calculate the optimal n and transform it back to the mag¬ 
netization m at the end. 

With increasing Dx/Dy, we see that the SkX and spi¬ 
ral phases become more stable relative to the vertical 
cone, and their region of stability extends into the easy- 
plane regime. To understand this, consider increasing the 
Rashba D±^ keeping D\\ fixed. The energy of cone m{z) 
depends only on Dy, and is unchanged as D± increases. 
In contrast, the SkX and spiral, with m = m(x,^), uti¬ 


lize the full D = (Dy + to lower their energy. 

Helicity and Ferrotoroidic Moment: We find 
that the spin textures smoothly evolve as a function of 
Dx/Dy. The spiral continuously changes from a Bloch- 
like (helical) spiral in the Dresselhaus limit to a Neel-like 
(cycloid) spiral in the Rashba limit. In between, the spins 
tumble around an axis at an angle P = tan“^(Dx/Dy) 
to the q-vector of the spiral. Similarly the skyrmions 
smoothly evolve from vortex-like textures in the Dressel¬ 
haus limit to hedgehogs in the Rashba limit, as seen in 
the insets of Fig. 1. In fact, 7 = 7 r /2 —;5is the “helicity”!^ 
of the skyrmions. 

Our results imply that D^/D\\ controls the helicity 7 , 
where Rashba could be tunable by electric field at an 
interface or by strain in a thin film. The ability to tune 7 
could be important in several ways. There is a recent pro¬ 
posal to use helicity to manipulate the Josephson effect 
in a superconductor/magnetic-skyrmion/superconductor 
junction^^. Another interesting phenomenon that de¬ 
pends on the helicity of skyrmions ^ is th e “ferrotoroidic 
moment” t = (1/2) j dPr[Y x m(r We will show 

elsewhere that t = to sinyz for the SkX. 

Rashba limit: Next we turn to the D\\ = 0 results 
in the right panel of Fig. 1, where one has the maximum 
regime of the stability for the spiral and the hexagonal 
SkX, in addition to a small region with a square lattice 
SkX (first predicted in ref. [25]) . an elliptic cone phase 
and a tilted F M. This phase diagram improves upon all 
previous work^l^^^^^ as explained in detail in the Sup¬ 
plementary Material. 

The tilted FM, which spontaneously breaks the U{1) 
symmetry of R (in a field), has ruz = Hj^A and exists in 
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the regime 2A> H and AJjD^ > 2 for D\\ = 0. We also 
see a new phase where the spins trace out a cone with 
an elliptic cross-section. The elliptic cone axis makes an 
angle Oq = /2A) with z, and the spatial variation 

of m is along a q- vector in the x-y plane. 

The nature of various phase transitions is discussed in 
the Supplementary Materials. In Fig. thick lines de¬ 
note continuous while thin lines denote first-order transi¬ 
tions. {A = 2, H = A)J/D^ is a Lifshitz poinfP^at which 
a state without broken symmetry (polarized FM) meets 
a broken symmetry (tilted FM) and a modulated (elliptic 
cone) phase. 

Let us next consider deviations from pure Rashba limit 
to see how the extreme right panel of Fig. 1 evolves into 
the 7 ^ 0 phase diagrams. As soon as one breaks bulk 
inversion, an infinitesimal leads to the tilted FM being 
overwhelmed by the vertical cone, which gains Dressel- 
haus DM energy. On the other hand, the elliptic and 
vertical cone states compete for Dy 7 ^ 0 and for some 
small, but finite, the vertical cone wins. 

Spin textures and topological charge: There are 
interesting differences between the skyrmions for H < 2A 
and H > 2A. {H = 2A is marked as a dashed line in the 
phase diagrams of Fig. 1 ). First, let us look at the spin 
textures. For SkX with H > 2 A, which have been the 
focus of all the past work, the spins at the boundary of 
the unit cell (u.c.) are all up, parallel to the field. Hence 
one can think of isolated skyrmions in a fully polarized 
FM background; see Fig. 2 left-most panels. It is the 
identification of the point at infinity in real space for an 
isolated skyrmion that lets us define a map from 
and use the homotopy group 7 r 2 ( 5 '^) = Z to characterize 
the topological charge or skyrmion number 

In contrast, when H < 2A, we find that the spins at the 
boundary are not all pointing up and the only constraint 
is periodic boundary conditions on the u.c.; see Fig. 2. 
There is no way to isolate this spin texture in a FM 
background. We must now consider the map r ^ ni(r) 
from the u.c., which is a 2 -torus to 5'^ in spin space. 
(Such maps are well known when represents a Bril- 
louin zone in k-space, but the mathematics is identical.) 
This map is characterized by an integer Chern number 
= /u.c. where x(r) = m • x dym)/4TT 

is the topological charge density. In fact, one can use this 
definition of A^sk for all values of iL/ 2 A. 

From the A = 0 panel on the left side of Fig. 2 , we 
see that x(r) is concentrated near the center of the u.c. 
and it is always of the same sign, as it is for all H > 2A. 
With increasing A, once H < 2 A, we see that x(r) begins 
to spread out and even changes sign within the u.c. In 
the square SkX phase x is again concentrated, but this 
time in regions near the center and the edges of each u.c. 
along with regions of opposite sign at the u.c. corners. 
For H ^ 2 A, the spin textures in the SkX phases are es¬ 
sentially composed of vortices and anti-vortices. Never¬ 
theless, the Chern number argument shows that the total 
topological charge in each u.c. is an integer; A^gk = “1 
in all the panels of Fig. 2. 


Discussion: Previous theories on understanding 

the increased stability on skyr mion s in thin films 
of non-centrosymmetric material d^^ * ^^^ ^ focussed pri¬ 
marily o n the ch anges in uniaxial magnetocrystalline 
anis otrop > h^ l ^^ l ^^ -^ with thickness, or on finite-size ef- 
fects^SI^. In fact, the latter can give rise to spin-textures 
more complicated than skyrmion crystals, with variations 
in all three directions. However, none of these theories 
take into account the role of broken surface inversion and 
Rashba SOC. As we have shown here, a non-zero Rashba 
D± leads to a greatly enhanced stability of the SkX 
phase, particularly for easy-plane anisotropy, while at the 
same time giving a handle on the helicity of skyrmions 
with interesting internal structure. 

We note that the phase diagrams in Fig. apply to all 
systems with broken mirror symmetry, with or without 
bulk inversion. Mirror symmetry can be broken by cer¬ 
tain crystal structures in bulk materials, by straiiP^ in 
thin films, or by electric fields at interfaces. For systems 
with D|| = 0 the vertical cone phase, which dominates 
much of the phase diagram for 7 ^ 0 , simply does not 
exist. After our paper was written, we became aware of 
the very recent observatiorP^ of hedgehog-like skyrmions 
in the magnetic semiconductor GaV 4 S 8 , a polar mate¬ 
rial with broken mirror symmetry that is dominated by 
Rashba SOC. Skyrmions are, however, stabilized in this 
material only at finite temperature due to the large easy- 
axis anisotropy. 

In conclusion, we have made a comprehensive study of 
the T = 0 phases, with a focus on skyrmion crystals in 
chiral magnets with two distinct DM terms. arises 
from Rashba SOC and broken surface inversion, while 
comes from Dresselhaus SOC and broken bulk inver¬ 
sion symmetry. We predict that increasing the Rashba 
SOC, via strain or electric field, and tuning magnetic 
anisotropy towards the easy-plane side will greatly help 
stabilize skyrmion phases in thin films, surfaces, and in¬ 
terface magnetism. Our results are very general, based 
on a continuum “Ginzburg-Landau” energy functional 
whose form is dictated by symmetry. We hope that it 
will motivate ab-initio density functional theory calcu¬ 
lations of the relevant phenomenological parameters en¬ 
tering our theory and an experimental investigations of 
skyrmions in Rashba systems. 
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Appendix A: Continuum Free Energy 

In this appendix we first introduce the continuum free 
energy functional that we use to model the spin textures 
in a chiral magnet. The free energy for a magnetic sys¬ 
tem with broken bulk inversion and mirror symmetries is 
F[m(r)] = f d^rJ^(m(r)) where 

J’(m(r)) = (J/2)(Vm)^ (Al) 

+D|| m • (V X m) 

-i-D± m • ((z X V) X m) 

-\-Aml — HrUz. 

and (Vm)^ is shorthand for ^-axis 

is the axis of broken mirror symmetry. Here J is the 
exchange stiffness, x is the DM 

vector, A is the magnetic anisotropy, and H is the field. 
We measure all the lengths in units of microscopic lattice 
spacing a, which we set to unity, so that J, Dy, A 
and iL all have units of energy. We are interested in low 
temperature behavior so we ignore fluctuations in the 
magnitude of the local magnetization m and impose the 
constraint that |m(r)| = 1. 

In our free energy functional we take the normal to the 
plane in which mirror symmetry is broken and the easy- 
axis direction to be the same, namely z. For simplicity, 
we also choose the external field to be along the same 
direction. One can, of course, imagine more general situ¬ 
ations in which these directions are not all the same, but 
the “simple” case treated here has sufficient complexity 
that it must be thoroughly investigated first. 

Next we define parameters D and /3 such that 

D\\ =Dcos/3 and Dj_ =Dsmp. (A2) 

It is convenient to rewrite using the natural energy 
and length scales in the problem. We measure energies 


in units of D‘^/J and lengths in units of JjD. In scaled 
variables the free energy density is given by 

J’(m(r)) = (l/2)(Vm)^ (A3) 

+m • ([cosyd V + sin/3z X V] x m) 
-f-Arnl — Hrriz^ 

which depends on three dimensionless variables AJjD^^ 
HJID‘^ and [3 = tan“^(Dx/D||). In the main paper, 
we explicitly write the anisotropy and field as AJjD^ 
and but in the Appendices we simplify notation 

and denote them as just A and H. (The total energy F 
depends on an inconsequential overall factor of {J/D)^ 
coming from the integration over the volume.) 

Next we briefly discuss the phases that we find as a 
function of A, H and /3. The two ferromagnetic (FM) 
phases - fully polarized and tilted - are states with no 
spatial variations. The (vertical) cone phase is a non- 
coplanar state which has only z-variations, along the 
magnetic field direction, so that m = m(^). These states 
can be treated analytically, as discussed in Appendix 

The spiral and skyrmion phases have a local magne¬ 
tization of the form m = m(x,^), as does the elliptic 
cone phase. Specifically, the spiral has spatial variation 
along a single direction, say x, in the plane perpendicu¬ 
lar to the field. The SkX phases have magnetic texture 
varying in both x and y. We use numerical methods for 
the analysis of phases with m = m(x, ^) as described in 
Appendix 0 Note that we do not consider states where 
m has non-trivial variations in all three coordinates; see, 
e.g., ref. [301 

Rotation: We next give the details of a transforma¬ 
tion (introduced in the main text) that greatly simplifies 
the analysis for states where m = m(x, y). We make the 
rotation 

( cos (3 sin yd 0 \ 

— siuyd cosyd 0 \ m = Rz{—l3)m. (A4) 

0 0 l) 

where (3 = tan~^(D x/-P ||)- Expressed in terms of m the 
free energy density ( |A3| ) simplifies to 

T = (l/2)(Vn)2 -1- n . (V X n) + (A5) 


provided that m, and thu s n, depends only on x and 
but not on 2 ;. This result (A5) has the same form as the 
free energy density in the pure Dresselhaus limit. After 
solving the problem in the n representation, we must 
transform back to m to find the actual spin texture. 

It is easy to see that the the exchange term is invari¬ 
ant under m ^ n, i.e., (Vm)^ —(Vn)^, as are the 
anisotropy term and Zeeman term coupling to H. The 
only term that has a nontrivial transformation is DM 
term in (|A3|). We symbolically write the DM term as as 
where 


V = cos yd V -f sin yd (z X V) 

( cos (3dx — sin f3dy \ 
sin pdx + cos /3dy j . ( A6) 

0 / 
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Here we set Vz = cos /3dz 0, because we focus on tex¬ 
tures that have no ^-variation, as already stated above. 
A straightforward calculation, using m = Rz{P) n, then 
allows us to derive 


m • (P X m) = n • (V X n) 


(A7) 


which in turn leads to eq. (A5). 


There is a slick way to obtain this same result by noting 
that the Free energy is invariant under a combined rota¬ 
tion in spin-space and in real space about the 2 ;-axis. (A 
combined spin and spatial rotation is needed because of 
SOC. and the 2 ;-axis is singled out by broken mirror sym¬ 
metry). In fact, it is simple to proceed with the general 
case where we reta in Vz = cos /3dz . The transformation 
Rz{—l3) of eq. (A4) acts only in spin-space. Thus to write 
the free energy in terms of n, we need to also rotate V 
in real space, so that V transforms to 

cos /3[cos I3dx — sin ^dy] + sin [sin ^dx + cos ^dy] 

— sin /3[cos pdx — sin ^dy] -h cos ;d[sin ^dx + cos ^dy] 
cos pdz 


This can be simply written as 


P ^ [V — (1 — cos P)dz ^]. (A8) 

Thus, for any magnetic texture m(x,^, z)^ the DM term 
can we written in general as 


+ n • (V X n) — (1 — cos /3)n • {dz{^ x n)) (A9) 


For a magnetic texture m(x, y) that does not vary along 
the z-axis the£-derivative terms vanish and this result 
simplifies to (A5) derived above. 


Appendix B: Ferromagnetic and Cone Phases 

In this appendix Q we consider phases which can 
be treated analytically: the polarized ferromagnet (FM), 
tilted FM and the vertical cone phase. 

Ferromagnets: The free energy density for a FM 
state is 


T = Arv?^ — HrUz. (Bl) 

Minimizing the free energy is easy in this case; we need 
to solve 


dT 

0 = -— = 2Amz - H (B2) 

oruz 

along with the constraint |mp = 1. The solution is 


m 


* 

2; 


r 1 H> 2 A 


^ J'* 


A-H H>2A 


(B3) 

We call the solution with m* = 1 a polarized FM since 
the magnetization is aligned with the magnetic field. The 


solution with m* < 1 we call tilted FM since the mag¬ 
netization is tilted away from the magnetic field. Unlike 
the polarized FM, the tilted FM spontaneously breaks 
the U{1) symmetry of spin rotation around 2 : axis in Bl 
Cone: Another simple class of magnetic states are 
textures 01 ( 2 :) that do not break translational symme¬ 
try in the x-y plane. We will find that the optimum 
configuration is a cone texture. This texture is called a 
cone because the magnetic moments trace out a cone as 
a function of 2 :; see the illustration in figure [BTj 

For textures m(^) with tra nslat ion symmetry in the 
x-y plane the Rashba term in (A3), with strength sin/3, 


does not contribute to the free energy. To implement the 
constraint |m( 2 :)| = 1 we define angular variables 0{z) 
and (j){z) such that 

m(^) = (cos (j) sin sin 0 sin 0, cos 0 ). (B4) 

In terms of 0{z) and ^( 2 :) the free energy density is 

R{0{z)^(j){z)) = ^(^0^ + 0 — (j)' cos P sin^ 0 


+A cos^ 0 — H cos 0. 
The Euler-Lagrange equation for ^( 2 :), 


(B5) 


f) f)T f) 

can be integrated with the result 

0' sin^ 0 — cos p sin^ 0 = C. (B7) 

Thus the free energy density can be written 

1 , 


H0{z)) = -{e'Y + f{e) 


(B8) 


where / does not depend on 6'. If is a minimum of / 
then f{0{z)) > f{0*) for all It is clear that 0{z) = 
is an extremum for the free energy obtained from (B5). 


Given that a constant 0{z) = minimizes the free en¬ 
ergy it is easy to find the optimum value 0* and the 
extremal function (j) which is 


^*( 2 :) = ^ cos j3. 


(B9) 


There are two solutions for . The first solution is ^* = 0 
which is a ferromagnetic solution with energy F = A — H. 
The second solution, with 


6>* = cos ^{H/{2A + cos^ (3)), 


(BIO) 


is called a cone texture. To make contact with the 
main text, we must recall that lengths are measured 
in units of J/D, so that eq. ( |B9| ) becomes 0*(^) ^ 
z{D/J){D\\/D) = zD \\/and energies A and H in units 
of DVJ, so that r ^ co^-\H/{2AfD\/J)]. 

We will occasionally refer to the cone as a vertical cone 
to distinguish it from the elliptic cone found in Appen¬ 
dices and The elliptic cone varies in the x-y plane 
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FIG. B.l. Illustration of the vertical cone phase, with 
q-vector (red arrow) along the z-axis, and the elliptic cone 
phase. In the elliptic cone phase the magnetization traces out 
an elliptic cone, i.e., the cross section is an ellipse rather than 
a circle. The elliptic cone phase shown here is for the = 0 
limit. If the cone height is decreased so that the spins lie in 
the x-z plane the conhguration becomes a cycloid (Neel-like 
spiral). 


while the vertical cone varies along the z-axis. An illus¬ 
tration of the vertical and elliptic cone textures is given 
in figure B.l In the Rashba limit, D\\ = 0, the cone 


phase is not stable for any values of H and A. For finite 
the tilted FM phase is not stable anywhere and the 
vertical cone takes its place in the phase diagram. 


Appendix C: Numerical Methods 

In appendix 0 we consider phases that cannot be 
treated analytically: the spiral, elliptic cone, square 
skyrmion crystal and hexagonal skyrmion crystal phases. 
For these phases the free energy density needs to be inte¬ 
grated and minimized numerically; we use conjugate gra¬ 
dient minimization to achieve this. All of these phases 
are of the form in(x,^) so we use the transformed free 
energy density ( |A5| ) and we will refer to n(x, y) as the tex¬ 
ture in the spin rotated frame. To implement numerical 
integration of the free energy some boundary conditions 
need to be specified. We use three boundary conditions: 
periodic along the x-axis, periodic with square symmetry 
and periodic with hexagonal symmetry. 

For each boundary condition the minimization proce¬ 
dure is very similar: convert the function n(x, y) to a 


vector s(i,jf), write the integral as a sum, then minimize 
to find the optimum vector. We discuss this procedure 
in detail only for the simplest case (periodic along the 
x-axis) and we discuss the important differences for the 
other two cases. 

Spiral and elliptic cone: First we consider textures 
which are periodic along a single axis; here we choose the 
X axis without loss of generality. The periodic boundary 
condition means that the texture n(x) satisfies n{x) = 
n{x + L). Since the free energy density is uniform along 
the y and z axes the free energy is just a one-dimensional 
integral along the x axis. 

To evaluate this integral numerically we can convert it 
to a sum. This involves discretizing the function n(x); 
let s(i) = n((i — l)Ax) where Ax = L/{N — 1). Then 
s(i +A) = s(i) is the periodic boundary condition. Next 
we replace all derivatives in the free energy density with 
finite differences, i.e., dxn{x) (s(i + 1) — s{i))/Ax for 
the point x = {i — l)Ax. Lastly we replace the integral 
with a sum (f^ dx Ax^f^-y). At this point the free 
energy can be computed given a configuration s{i) and 
this will be a good approximation to the true free energy 
when N is large, i.e., F(s) F[n] is a good approxima¬ 
tion when N is large. 

We also have the constraint \n{x)\ = 1. In terms 
of s the constraint is |s(i)| = 1. We impose the con¬ 
straint by introducing angular variables, 0{i) and (j){i)^ 
at each site so that s{i) = (cos 0 sin sin 0 sin cos ^). 
Now we can think of F as a function of the 2N compo¬ 
nent vector {6>(1), 0{N), •••5 0(^)} = {^5 0} = P- 

We can minimize the vector function F(p) using con¬ 
jugate gradient minimization. Conjugate gradient mini¬ 
mization accepts as its input a function, in this case F, 
the gradient of that function VpF, and an initial vector, 
call it po; the output of conjugate gradient minimiza¬ 
tion is a local minimum p*. We choose single-q spirals 
m(x) = (cosgx, singx, 0) for our initial condition where 
q = 21: jL. This choice of initial condition does not limit 
the output of conjugate gradient minimization to single-q 
spirals; the output has a more complicated Fourier space 
structure in general. 

Note that the function F(p) must be a real valued 
function of p. In particular we need to specify A, H and 
the periodicity L for F to be a real valued function. Note 
that L is also a variational parameter in this case. In 
order to find the optimum value of L we include L in the 
conjugate gradient minimization procedure, i.e., we think 
of F(p, L) as a function of the 2N + 1 component vector 
{p,F}. We then sweep across the phase diagram and 
find the optimum configuration for each value of —1.6 < 
A < 3.1 and 0 < F < 4.1 in steps of 0.1 along both axes. 

The results of this minimization procedure are incor¬ 
porated into the phase diagram in Fig. 1 of the main 
text. The phases that we find are spiral phases and el¬ 
liptic cone phases. Since we minimize the simplified free 
energy ( |A5[ ) the spirals we find are all Bloch-like spi¬ 
rals. These textures can be expressed in terms of a single 
function 0{x) where n{x) = (0, sin^(x), cos^(x)). The 













textures obtained in the lab frame are related to these 
texture by a rotation by P about the z-axis as discussed 
in Appendix[A| i.e., m(x) = (sin/3 sin cos/3 sin cos ^). 

The other phase we find is the elliptic cone. This phase 
is more complicated than the spiral phase. Where the spi¬ 
ral can be described in terms of a single function, 6 >(x), 
the elliptic cone requires two functions in general, i.e., 
0{x) and (j){x). As a consequence of the (j){x) dependence 
all components of the magnetic moment vary as a func¬ 
tion of X, not just the y and 2 ;-components like the spiral. 
An illustration of the elliptic cone is given in Fig. |B.l| and 
in Appendix we discuss a variational state that cap¬ 
tures the physics of the elliptic cone. 

Square skyrmion crystal: Next we consider tex¬ 
tures that are periodic along both axes, i.e., n(x,^) = 
n(x + 1/,^) = n(x,^ + L). We further restrict our¬ 
selves to textures with C 4 symmetry and we use this 
symmetry to reduce the number of points in the unit 
cell by a factor of 4 (see Fig. C.l for an illustration). 
We convert the free energy integral over the square unit 
cell to a sum and we convert derivatives to finite dif¬ 
ferences. Here we need a 2N x N component vector 

to describe a given magnetic texture n(x,^). The unit 
cell has x N sites. We need to specify the function, 
F, its gradient, VpF, and an initial vector po. We use a 
single-q vortex-like skyrmion as our initial configuration. 


I.e., 


n(x, y) = [-{y/r)± -h {x/r)y] sin(gr) + z cos(gr), (Cl) 

where r = x‘^ -\- y‘^ and we use q = tt jL. Note that this 

configuration is actually a square skyrmion crystal since 
we are applying square periodic boundary conditions. As 
we have already stressed, the simple form of the initial 
configuration does not limit the form of the output of 
conjugate gradient minimization. 

To eliminate discretization error and approach the con¬ 
tinuum limit, we extrapolate our results to the N ^ 00 
limit to obtain the phase diagram in Fig. 1 of the main 
text. We use a polynomial fit of the energy vs Ax for 
N = 40,..., 60 and we find that the Ax = 0 limit is the 
same whether we fit to a polynomial of degree 3, 4 or 5. 
This gives us confidence that we are sufficiently close to 
the N ^ oc limit that extrapolation is valid. 

The optimum configurations we find using conjugate 
gradient minimization are square skyrmion crystals. The 
topological charge density n • {dxn x dyu) is not con¬ 
centrated at the center of each unit cell like it is for 
skyrmions found in the Dresselhaus limit. The topo¬ 
logical charge density is concentrated at four points in 
the unit cell: the center, corners and edges of the unit 
cell (see Fig. 2 in the main text) with opposite sign at 
the corners. Such configurations have been referred to 
as meron crystal^^ nevertheless, the skyrmion number 
for this phase is quantized in integer units by the Chern 
number argument given in the main text. The square 
skyrmion crystal phase is stable in a small pocket of the 
F|| = 0 phase diagram. 



FIG. C.l. Illustration of square (left) and hexagonal (right) 
unit cells. The dark lines indicate the region of independent 
spins (indicated by red dots). Spins in the dashed region (gray 
dots) are determined by C 4 symmetry for the square system 
and C3 symmetry for the hexagonal system. 


Hexagonal skyrmion crystal Lastly we consider 
textures with hexagonal symmetry, i.e., x(x,^) = m{x-\- 
ai,^ + a 2 ) = m(x + 6 i,^ + 62 ) where (ai,a 2 ) and ( 61 , 62 ) 
are lattice vectors for a regular triangular lattice with lat¬ 
tice spacing L and we also assume C 3 symmetry. In figure 
|C.l| we illustrate how we use the Cs symmetry to reduce 
the number of points in our simulation by a factor of 3. 
Besides the symmetry of the problem, the minimization 
follows exactly as for the square symmetry case. 

The optimum configurations we find using conjugate 
gradient minimization are hexagonal skyrmion crystals. 
The hexagonal crystal has a much large range of stabil¬ 
ity than the square crystal. Above the line H = 2A the 
hexagonal crystal has topological charge density concen¬ 
trated in the center of each unit cell and the topological 
charge density is essentially zero throughout the unit cell. 
All skyrmions found in the Dresselhaus limit are of this 
type. Below the H = 2A line the topological charge den¬ 
sity becomes less dense at the center of the unit cell and 
the topological charge is smeared across the unit cell. 
Close to the square crystal phase the topological charge 
density gathers at the center, corners and edges of the 
unit cell, with the skyrmion number still remaining quan¬ 
tized. 


Appendix D: Variational Solution 

In appendix we discuss a variational state that can 
be integrated exactly and captures the physics of the po¬ 
larized FM, tilted FM, spiral, elliptic cone and vertical 
cone phases. The main point of this variational state is to 
gain analytical insight into the new elliptic cone phase; 
in particular, this variational state confirms the phase 
boundary between the elliptic cone and the tilted FM is 
A = 2 and the phase boundary between the elliptic cone 
and the polarized FM is iL = 2A. 

Consider the spin-texture 

m(r) = eimi cosq • r + e 2 m 2 sinq • r + esms{T) (Dl) 
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where ei, 62 an d 63 are orthonormal vectors. The func¬ 
tion m3(r) = — (mi cosq • r)^ — (m2 sinq • r)^ is de¬ 

fined to enforce the normalization |m(r)| = 1 at ev¬ 
ery point. We must restrict mi and m 2 to the range 
— 1 < mi, m2 < 1 to make sure m3 is a real number. No¬ 
tice that this state varies along the q axis and is uniform 
perpendicular to this axis. 

Limits: Here we discuss various limits of the varia¬ 
tional state. The simplest case is the limit mi = m2 = 0 
where the variational state reduces to m = 63 which 
is just a ferromagnetic configuration, either polarized 
(e 3 = z) or tilted. 

Another simple limit is mi = m2 = 1. In this case 
the variational state becomes m = ei cos q • r + 62 sin q • r 
which is a spiral state. Notice that this is a very restricted 
state which is a delta function in Fourier space. In gen¬ 
eral spiral textures have complicated Fourier structure so 
we do not expect this state to capture the spiral phase 
boundaries quantitatively, only qualitatively. There is a 
special point in the phase diagram [H = 0 and A = 0) 
where the spiral actually has this simple form. 

When 0 < mi = m2 < 1 the state is a circular cone, 
i.e., the magnetic moments trace out a cone in spin space 
as a function of ^ when q^ = qz. The circular cone found 
in Appendix [B] is exactly of this form with q = cos/3z, 
mi = sin6>*, ei = x and 62 = z. 

Finally we discuss the most complicated limiting case 
of thi s va riational state. When 0 < mi < m2 < 1 the 
state ( |D1| ) is an elliptic cone, i.e., the magnetic moments 
trace out an elliptic cone as a function of x, for example, 
when q = X. An elliptic cone has a base and a tip, similar 
to a cone, but the cross sections are ellipses rather than 
circles. The actual elliptic cone found using numerics is 
not as simple as the variational state. In particular the 
numerical solution does not have simple Fourier structure 
along the ei and 62 axes, in contrast to the variational 
state |Dl[ However, near the polarized FM-elliptic cone 
and tilted FM-elliptic cone phase boundaries the elliptic 
cone becomes more and more like the variational state 
|dTJ so the phase boundary found using this variational 
state is quantitatively correct. 

Exact integration: It is convenient to trade the 
for Euler angles. This can be achieved by writing 


ei = nMny{02)nM±i (D2) 


where xi = x, X 2 = y and X 3 = z. This state ( |D1| ) 
has 8 variational parameters: 6>i, 6^2, 6^3, mi, m2, and 
q. Without loss of generality we can choose q to lie in 
the x-z plane leaving 7 variational parameters. The free 


energy in terms of these 7 parameters is 
F = Fj Fd + Fa + Fh 
Pj = + il)iX - -”^ 2 )) 

Fd = —mim2 cos I3{qz cos O 2 + qx cos 0^ sin O 2 ) 
—mi m2 sin Pqx sin O 2 sin 0^ 

Fa= y cos^ 6*2(2 - TOi -TO2) 

+ — sin^ 6>2((mi cos6>i)^ + (m 2 sin6>i)^) 


(D3) 


Fh = — 


H cos O 2 
27r 


p27T 

/ — (mi cosi^)^ — (m2 sini^)^. 

Jo 


Fh can be expressed in terms of the complete elliptic 
integral E{x) = \/l — xsin^ OdO, i.e.. 


H 

1^12 


Fh = —^ a/I ~ ^2^ ( — 


(mi - m2)(mi + m 2 ) 

I - ml 


The free energy has been integrated analytically but the 
variational parameter minimization involves transcen¬ 
dental equations. Numerically solving these equations 


produces the phase diagram in Figure D.l in the = 0 
limit. 

Important features of this phase diagram are the po¬ 
larized FM-elliptic cone and tilted FM-elliptic cone phase 
boundaries. These boundaries agree to arbitrary preci¬ 
sion with the phase boundaries found using numerics, 
i.e., H = 2A for the polarized FM-elliptic cone transition 
and A = 2 for the tilted FM-elliptic cone transition. The 
spiral phase boundary is only qualitatively correct due to 
the simple Fourier structure of the variational state. 


Appendix E: Phase Transitions 


In this Appendix we discuss the nature of the various 
phase transitions shown in Fig. 1 of the main text, where 
thick lines indicate continuous transitions and thin lines 
indicate first order transitions. All phase transitions are 
determined by com parin g energy curves near the phase 
boundary (see Fig. E.l). Energy curves that have the 
same slope on both sides of the phase boundary cor¬ 
respond to continuous transitions, discontinuous slopes 
correspond to first order transitions. In the = 0 limit 
of the phase diagram we find only one continuous phase 
transition, between the polarized EM and the cone (ver¬ 
tical). In the = 0 limit we find four continuous phase 
transitions, three of which meet at a point called a Lif- 
shitz poinfP^. Below we give reasoning for the continuous 
or first-order nature of the various transitions. 

Dresselhaus limit: The D± = 0 limit of the phase 
diagram has been previously studiecP^. The phase tran¬ 
sitions are all first order except the polarized EM to cone 
transition, which is continuous. At this transition the 
cone radius goes continuously to zero. 
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FIG . D.l. Phase diagram obtained from variational Ansatz 
(Dl) for = 0. Thick lines denote continuous transitions 
and thin lines denote first order transitions. The variational 
state does not allow for skyrmion phases so we do not expect 
them in this phase diagram. The main result is that the phase 
boundaries at = 2A and AJjD^ = 2 agree to arbitrary 
precision with the numerical results, i.e., the polarized FM- 
elliptic cone and tilted FM-elliptic cone phase boundaries are 
identical to the boundaries in the T)|| = 0 limit of Fig. 1 of 
the main text. 


HJ/D2=0.8 



HJ/D2=0.8 



The spiral-cone, spiral-skyrmion crystal and cone- 
skyrmion crystal transitions are all first order. In each of 
these transitions the phases on either side of the phase 
boundary have distinct broken symmetries. Generically 
we expect to find first order transitions between two 
phases with different broken symmetries. 

For the first order polarized FM-spiral phase transition 
the polarized FM phase has no broken symmetry while 
the spiral phase breaks translational symmetry. Further¬ 
more, the wavelength of the spiral remains finite at the 
phase boundary. 

Our numerical results for the polarized FM-skyrmion 
crystal show that the transition is first order; however, 
it is an unusual first order transition with a diverging 
length scale associated with the optimal SkX unit cell 

siz5^. 

Rashba limit: We discuss here only the phase transi¬ 
tions that are present in the = 0 phase diagram and 
are not present in the D± = 0 phase diagram. 

The elliptic cone transitions continuously to the polar¬ 
ized FM and tilted FM phases. Near the phase boundary 
the radius of the cone going continuously to zero. There 
is a point, A = 2 and H = 4, where these three phases 
meet. This is a Lifshitz poinfP^ at which a “symmetric” 
phase, in our case the polarized FM, meets a broken sym¬ 
metry phase, the tilted FM, and a spatially modulated 
phase, the elliptic cone; see Sec. 4.6 of ref. [26l 

The elliptic cone-spiral phase transition is also con¬ 
tinuous, in contrast with the vertical cone-spiral phase 
boundary. Near the elliptic cone-spiral phase boundary 


FIG. E.l. Free energy relative to the tilted FM (top) as 
a function of AJ/D^ at fixed HJ/D^ = 0.8 and the deriva¬ 
tive dFjdA of the free energy (bottom) for the hexagonal 
skyrmion crystal (red), square skyrmion crystal (blue), ellip¬ 
tic cone (black) and tilted FM (green) phases. Phase tran¬ 
sitions are marked by cyan lines. In the plot of dF/dA it is 
easy to see which phase transitions are continuous (elliptic 
cone-tilted FM) and which are first order (hexagonal-square 
skyrmion crystal and square skyrmion crystal-elliptic cone) 
by examining the jump discontinuities in the derivative of the 
free energy. The stable phase is indicated by a darker line in 
the bottom figure. 

the height of the cone goes continuously to zero. 

The elliptic cone-hexagonal skyrmion crystal, elliptic- 
cone square skyrmion crystal and hexagonal-square 
skyrmion crystal phase boundaries are all first order with 
distinct broken symmetries on each side of the phase 
boundaries. 


Appendix F: Rashba Limit Phase Diagram 

In appendixj^we comment in detail on how the Rashba 
limit phase diagram shown in the right-most panel of 
Fig. 1 of the main paper, goes beyond all previous works, 
as mentioned in the text. 

The H = ^ results for A > 0 can be compared by 
taking the T = 0 limit of the finite temperature analysis 
of ref. |20l For = 0 and A > 0 they find three phases: 
a spiral, a cone, and an in-plane FM. These are same 
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as the H = 0 limit of our results with our tilted elliptic 
cone simplifying to a cone whose axis is horizontal. Our 
phase boundaries are more accurate, however, because, 
as emphasized in the text, we do not restrict attention 
to a single-Q spiral as in ref. [20l Further, these authors 
did not analyze H ^ 0. 

For 7 ^ 0, our own earlier worlP^ presented clear 
evidence for the development of nontrivial spatial struc¬ 
ture in the topological charge density, but in that paper 
we worked with a SkX variational ansatz that forced the 
spins on the unit cell boundary to be pointing along the 
field direction. Thus we did not have the variational free¬ 
dom to see skyrmions where we must use a Chern number 
to understand the quantization of topological charge. We 
see here the larger hexagonal SkX region than in ref. |T9| 
and also the square SkX. 

Recently, ref. [25] predicted a small region of stability 
for the square SkX in the A > 0 regime, and our results 
are consistent with theirs as far as this feature of the 
phase diagram is concerned. However, both refs. [TOland 
[25] missed the elliptic cone phase. 


interactions that lead to shape anisotropy. 

Here we comment on the anisotropy contribution of 
the SOC that leads to the DM terms in the free energy 
functional, following refs. [18] and (TQ] Rashba SOC, which 
gives rise to a DM term D± ^ Agoc linear in the SOC cou¬ 
pling constant Agoc, also gives rise to a compass-Kitaev 
anisotropy of the form 


— 4 , f 8'^ 8'^ -i- qy qy \ 


(Gl) 


This term is often ignored in the literature because 
^ small. However, this argument is flawed 

because the energetic contribution of the DM term is of 
order D\/J and hence of the same order as that of A^. 
In fact, one can show that for a large class of exchange 
mechanisms A±J/D\ 1/2 in the limit of weak SOC; 
see ref. [T9] and references therein. In the continuum limit, 
and ignoring higher order derivative terms, eq. (Gl) leads 
to This is an important easy-plane contribu¬ 

tion to the total anisotropy. 

What about the analogous term arising from Dressel- 
haus SOC? This is of the form 


Appendix G: Magnetic Anisotropy 

In appendix [G] we discuss how Rashba spin-orbit cou¬ 
pling gives rise to an easy-plane magnetic anisotropy. In 
the paper we have a magnetic anisotropy term J^a = 
Ami energy functional, where A is treated as 

a phenomenological parameter that can be either A < 0 
(easy-axis) or A > 0 (easy-plane). As noted there, 
many mechanisms contribute to A, including atomic SOG 
which gives rise to single-ion anisotropy and and dipolar 


- ^11 J (G2) 

with A||J/Dy =1/2. If we take the continuum limit, 
retain just the order terms and ignore higher order 
derivative terms, we get A||(m^ + + ml) which is 

just an additive constant of no consequence, since = 
1. So it is only the Rashba SOG that gives rise to the 
interesting easy-plane anisotropy, as one can see from 
symmetry alone. 
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